{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业内容\n",
    "\n",
    "> 基于阿尔法衰变半衰期，推导最小二乘法公式并用Latex在Markdown cell中打出来；\n",
    "\n",
    "> 导入数据，运用最小二乘法，使用前19个点拟合，并预测第20个点的值，与真实值进行比较。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 参考答案"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 1， 理论部分"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由公式：\n",
    "\n",
    "$$ \\log_{10} W = C - \\frac{D}{\\sqrt{E_\\alpha}} $$\n",
    "\n",
    "定义：\n",
    "\n",
    "$$ y = \\log_{10}{W} $$\n",
    "\n",
    "$$ x = 1/\\sqrt{E_{\\alpha}} $$\n",
    "\n",
    "记:\n",
    "\n",
    "$$C = b, -D = a$$\n",
    "\n",
    "所以有：\n",
    "\n",
    "$$ y = ax + b $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Case 1**\n",
    "\n",
    "使用L2 Loss (MSE)，有：\n",
    "\n",
    "$$ l(a, b) = \\sum_{i}^{n}{(y_i - \\hat{y_i}})^2 $$\n",
    "\n",
    "其中，\n",
    "\n",
    "$$ \\hat{y_i} = ax_i + b $$\n",
    "\n",
    "分别求损失函数对两个参数的偏导数，有：\n",
    "\n",
    "$$ \\Lambda_a = \\frac{\\partial l}{\\partial a} = \\sum_{i}^{n}{2(y_i - ax_i - b)(-x_i)} $$\n",
    "\n",
    "$$ \\Lambda_b = \\frac{\\partial l}{\\partial b} = \\sum_{i}^{n}{2(y_i - ax_i - b)(-1)} $$\n",
    "\n",
    "令它们等于0，可得：\n",
    "\n",
    "$$ \\sum_{i}^{n}{y_i x_i} - a\\sum_{i}^{n}{(x_i)^2} - b\\sum_{i}^{n}{x_i} = 0 \\ \\ \\ \\ (1)$$\n",
    "\n",
    "$$ \\sum_{i}^{n}{y_i} - a\\sum_{i}^{n}{x_i} - b\\sum_{i}^{n}{1} = 0 \\ \\ \\ \\ (2)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对(1)式乘上$\\sum_{i}^{n}{1}$， 对(2)式乘上$\\sum_{i}^{n}{x_i}$，可得：\n",
    "\n",
    "$$ n\\sum_{i}^{n}{y_i x_i} - an\\sum_{i}^{n}{(x_i)^2} - bn\\sum_{i}^{n}{x_i} = 0 \\ \\ \\ \\ (3)$$\n",
    "\n",
    "$$ \\sum_{i}^{n}{y_i}\\sum_{i}^{n}{x_i} - a(\\sum_{i}^{n}{x_i})^2 - bn\\sum_{i}^{n}{x_i} = 0 \\ \\ \\ \\ (4)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "很显然地，有：\n",
    "\n",
    "$$ a = \\frac{n\\sum_{i}^{n}{y_i x_i}-\\sum_{i}^{n}{y_i}\\sum_{i}^{n}{x_i}}{n\\sum_{i}^{n}{x_i^2} - (\\sum_{i}^{n}{x_i})^2} $$\n",
    "\n",
    "$$ b = \\frac{\\sum_{i}^{n}{y_i} - a\\sum_{i}^{n}{x_i}}{n} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "特别地，注意到这里的式子里面有n，所以对a而言，分子分母同时除以一个n或者除以n平方，都可以得到平均值的表达式；对b而言，裂项之后的结果两项都是平均值。\n",
    "\n",
    "所以同样可以得到另外一套结果：\n",
    "\n",
    "$$ a = \\frac{\\bar{y_i x_i} - \\bar{y_i}\\bar{x_i}}{\\bar{x_i^2} - \\bar{x_i}^2} $$\n",
    "\n",
    "$$ b = \\bar{y_i} - a\\bar{x_i} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Case 2**\n",
    "\n",
    "以上是选用MSE的结果，如果选择MAE，则是：\n",
    "\n",
    "$$ l(a, b) = \\sum_{i}^{n}{\\sqrt{(y_i - \\hat{y_i}})^2} $$\n",
    "\n",
    "其中，\n",
    "\n",
    "$$ \\hat{y_i} = ax_i + b $$\n",
    "\n",
    "分别求损失函数对两个参数的偏导数，有：\n",
    "\n",
    "$$ \\frac{\\partial l}{\\partial a} = \\sum_{i}^{n}{\\frac{2(y_i - ax_i - b)(-x_i)}{|y_i - ax_i - b|}} $$\n",
    "\n",
    "$$ \\frac{\\partial l}{\\partial b} = \\sum_{i}^{n}{\\frac{2(y_i - ax_i - b)(-1)}{|y_i - ax_i - b|}} $$\n",
    "\n",
    "记：\n",
    "\n",
    "$$ |y_i - ax_i - b| = \\Lambda_i $$\n",
    "\n",
    "令它们等于0，可得：\n",
    "\n",
    "$$ \\sum_{i}^{n}{y_i x_i / \\Lambda_i} - a\\sum_{i}^{n}{x_i^2 / \\Lambda_i} - b\\sum_{i}^{n}{x_i / \\Lambda_i} = 0 \\ \\ \\ \\ (5)$$\n",
    "\n",
    "$$ \\sum_{i}^{n}{y_i / \\Lambda_i} - a\\sum_{i}^{n}{x_i / \\Lambda_i} - b\\sum_{i}^{n}{1 / \\Lambda_i} = 0 \\ \\ \\ \\ (6)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "同样地，对(5)式乘上$\\sum_{i}^{n}{1 / \\Lambda_i}$， 对(6)式乘上$\\sum_{i}^{n}{x_i / \\Lambda_i}$，可得：\n",
    "\n",
    "$$ \\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{y_i x_i / \\Lambda_i} - a\\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{x_i^2 / \\Lambda_i} - b\\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{x_i / \\Lambda_i} = 0 \\ \\ \\ \\ (7)$$\n",
    "\n",
    "$$ \\sum_{i}^{n}{x_i / \\Lambda_i}\\sum_{i}^{n}{y_i / \\Lambda_i} - a(\\sum_{i}^{n}{x_i / \\Lambda_i})^2 - b\\sum_{i}^{n}{x_i / \\Lambda_i}\\sum_{i}^{n}{1 / \\Lambda_i} = 0 \\ \\ \\ \\ (8)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "显然，可得：\n",
    "\n",
    "$$ a = \\frac{\\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{y_i x_i / \\Lambda_i} - \\sum_{i}^{n}{x_i / \\Lambda_i}\\sum_{i}^{n}{y_i / \\Lambda_i}}{\\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{x_i^2 / \\Lambda_i} - (\\sum_{i}^{n}{x_i / \\Lambda_i})^2} $$\n",
    "\n",
    "$$ b = \\frac{\\sum_{i}^{n}{y_i / \\Lambda_i} - a\\sum_{i}^{n}{x_i / \\Lambda_i}}{\\sum_{i}^{n}{1 / \\Lambda_i}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 2，编程部分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nucleus</th>\n",
       "      <th>Ealpha</th>\n",
       "      <th>tau_half</th>\n",
       "      <th>tau_units</th>\n",
       "      <th>tau_second</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Po206</td>\n",
       "      <td>5.22</td>\n",
       "      <td>8.80</td>\n",
       "      <td>day</td>\n",
       "      <td>7.603200e+05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Po208</td>\n",
       "      <td>5.11</td>\n",
       "      <td>2.90</td>\n",
       "      <td>year</td>\n",
       "      <td>9.145440e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Po210</td>\n",
       "      <td>5.31</td>\n",
       "      <td>138.00</td>\n",
       "      <td>day</td>\n",
       "      <td>1.192320e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Po212</td>\n",
       "      <td>8.78</td>\n",
       "      <td>0.30</td>\n",
       "      <td>microsecond</td>\n",
       "      <td>3.000000e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Po214</td>\n",
       "      <td>7.68</td>\n",
       "      <td>164.00</td>\n",
       "      <td>microsecond</td>\n",
       "      <td>1.640000e-04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Po216</td>\n",
       "      <td>6.78</td>\n",
       "      <td>0.15</td>\n",
       "      <td>second</td>\n",
       "      <td>1.500000e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>U228</td>\n",
       "      <td>6.69</td>\n",
       "      <td>9.10</td>\n",
       "      <td>month</td>\n",
       "      <td>2.358720e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>U230</td>\n",
       "      <td>5.89</td>\n",
       "      <td>20.80</td>\n",
       "      <td>day</td>\n",
       "      <td>1.797120e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>U232</td>\n",
       "      <td>5.32</td>\n",
       "      <td>72.00</td>\n",
       "      <td>year</td>\n",
       "      <td>2.270592e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>U234</td>\n",
       "      <td>4.77</td>\n",
       "      <td>247000.00</td>\n",
       "      <td>year</td>\n",
       "      <td>7.789392e+12</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  nucleus  Ealpha   tau_half    tau_units    tau_second\n",
       "0   Po206    5.22       8.80          day  7.603200e+05\n",
       "1   Po208    5.11       2.90         year  9.145440e+07\n",
       "2   Po210    5.31     138.00          day  1.192320e+07\n",
       "3   Po212    8.78       0.30  microsecond  3.000000e-07\n",
       "4   Po214    7.68     164.00  microsecond  1.640000e-04\n",
       "5   Po216    6.78       0.15       second  1.500000e-01\n",
       "6    U228    6.69       9.10        month  2.358720e+07\n",
       "7    U230    5.89      20.80          day  1.797120e+06\n",
       "8    U232    5.32      72.00         year  2.270592e+09\n",
       "9    U234    4.77  247000.00         year  7.789392e+12"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#load data\n",
    "data = pd.read_csv('./data/alpha_decay.csv')\n",
    "\n",
    "#unit conversion\n",
    "def UnitConversion(unit):\n",
    "    if unit not in ['microsecond', 'second', 'day', 'month', 'year']:\n",
    "        raise Exception('Invalid unit.')\n",
    "    conv_rate = 1\n",
    "    if unit == 'microsecond':\n",
    "        conv_rate = 1e-6\n",
    "    if unit == 'day':\n",
    "        conv_rate = 3600 * 24\n",
    "    if unit == 'month':\n",
    "        return 3600 * 24 * 30\n",
    "    if unit == 'year':\n",
    "        conv_rate = 3600 * 24 * 365\n",
    "    return conv_rate\n",
    "\n",
    "data['tau_second'] = [UnitConversion(item) for item in data['tau_units']] * data['tau_half']\n",
    "\n",
    "#check\n",
    "data.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x: [0.43768811 0.44237396 0.43396304 0.33748365 0.36084392 0.38404769\n",
      " 0.38662234 0.41204282 0.43355498 0.45786855 0.47192918 0.48795004\n",
      " 0.42640143 0.4397995  0.45175395 0.46324105 0.39872611 0.40422604\n",
      " 0.4152274  0.43073049]\n",
      "y: [ -6.04017095  -8.12037914  -7.23556737   6.36370421   3.62598161\n",
      "   0.6647342   -7.53185093  -6.41375162  -9.51531364 -13.0506781\n",
      " -15.03637905 -17.31215769  -9.5924796  -11.47620704 -13.23662036\n",
      " -15.56107113  -6.52382308  -7.30787589  -8.90349381 -11.39834384]\n"
     ]
    }
   ],
   "source": [
    "#define variables\n",
    "x = np.sqrt(1 / data['Ealpha'])\n",
    "W = np.log(2) / data['tau_second']\n",
    "y = np.log10(W)\n",
    "print('x: {}\\ny: {}'.format(x.values, y.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x2d340ca6978>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHSCAYAAAATyJnbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAWUElEQVR4nO3df4jt+X3X8df77trIlBQTsjFpkrmzkQTZWE3IEA0SUYkmDbQxjcWUQ20tMg20+If+YctVEMtFLEKptFCnUFE4MUggNv1hYgJNg2Bo7zWbuLGubjY7N5stuEmFFqemtPvxj3Mmd+7duT/mzpw5533O4wHDOd/P+fX93LMnz3y/853vqTFGAIBeLi17BQCA0xNwAGhIwAGgIQEHgIYEHAAaEnAAaOjhZa/Aca94xSvGzs7OslcDAC7M9evXvzbGeOS0j1upgO/s7OTatWvLXg0AuDBVdfAgj7MLHQAaEnAAaEjAAaAhAQeAhgQcABoScABoSMABoCEBB4CGBBwAGhJwAGhIwAGgIQEHgIYEHAAaEnAAaEjAAaAhAb+D6TTZ2UkuXZpdTqfLXiMAuOnhZa/AKppOk7295PBwtnxwMFtOkslkeesFAEdsgZ/gypWb8T5yeDgbB4BVIOAnuHHjdOMAcNEE/ATb26cbB4CLJuAnuHo12dq6dWxrazYOAKtAwE8wmST7+8nly0nV7HJ/3wFsAKwOR6HfwWQi2ACsLlvgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwv/OtGqeibJ7yX5oyR/OMbYXfRrAsC6u6jvA/8rY4yvXdBrAcDaswsdABq6iICPJP+pqq5X1d4FvB4ArL2L2IX+F8cYz1XVK5N8sqr+xxjjM0c3zqO+lyTb29sXsDoA0N/Ct8DHGM/NL/93ko8medttt++PMXbHGLuPPPLIolcHANbCQgNeVd9aVS89up7kryd5YpGvCQCbYNG70P9kko9W1dFrfWiM8fEFvyYArL2FBnyM8XSSP7fI1wCATeTPyACgIQEHgIYEHAAaEnAAaEjAAaAhAQeAhgQcABoScABoSMABoCEBB4CGBBwAGhJwAGhIwAGgIQEHgIYEfMNNp8nOTnLp0uxyOl32GgFwPxb6feCstuk02dtLDg9nywcHs+UkmUyWt14A3Jst8A125crNeB85PJyNA7DaBHyD3bhxunEAVoeAb7Dt7dONA7A6BHyDXb2abG3dOra1NRsHYLUJ+AabTJL9/eTy5aRqdrm/7wA2gA4chb7hJhPBBujIFjgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAsxGm02RnJ7l0aXY5nS57jQDO5uFlrwAs2nSa7O0lh4ez5YOD2XKSTCbLWy+As7AFztq7cuVmvI8cHs7GAboScNbejRunGwfoQMBZe9vbpxsH6EDAWXtXryZbW7eObW3NxgG6EnDW3mSS7O8nly8nVbPL/X0HsAG9OQqdjTCZCDawXmyBA0BDAg4ADS084FX17qp6sqqeqqofW/TrAcAmWGjAq+qhJD+b5DuTPJbk+6rqsUW+JgBsgkVvgb8tyVNjjKfHGH+Q5MNJ3rvg1wSAtbfogL8myVeOLT87HwMAzmDRAa8TxsYtd6jaq6prVXXt+eefX/DqAMB6WHTAn03yumPLr03y3PE7jDH2xxi7Y4zdRx55ZMGrAwDrYdEB/80kb6iqR6vqW5J8IMnHFvyaALD2FnomtjHGH1bVjyb5RJKHkvzCGOOLi3xNANgECz+V6hjjV5P86qJfBwA2iTOxAUBDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADS0sIBX1T+pqq9W1ePzn/cs6rUAYNM8vODn/6kxxr9Y8GsAwMaxCx0AGlp0wH+0qr5QVb9QVS876Q5VtVdV16rq2vPPP7/g1QGA9VBjjAd/cNWnkrzqhJuuJPlskq8lGUl+Ismrxxg/dLfn293dHdeuXXvg9QGAbqrq+hhj97SPO9PvwMcY77yf+1XVzyf55bO8FgBw0yKPQn/1scX3JXliUa8FAJtmkUeh/2RVvTmzXejPJPnhBb4WAGyUhQV8jPH9i3puANh0/owMABoScABoSMABoCEBB4CGBBwAGhJwAGhIwAGgIQEHgIYEHAAaEnBg40ynyc5OcunS7HI6XfYawekt8lzoACtnOk329pLDw9nywcFsOUkmk+WtF5yWLXBgo1y5cjPeRw4PZ+PQiYADG+XGjdONw6oScGCjbG+fbhxWlYADG+Xq1WRr69axra3ZOHQi4LChNvVI7Mkk2d9PLl9OqmaX+/sOYKMfR6HDBtr0I7Enk82YJ+vNFjgswbK3fh2JDf3ZAocLtgpbv47Ehv5sgcMFW4WtX0diQ38CDhdsFbZ+HYkN/Qk4XLBV2Pp1JDb0J+BwwVZl63cySZ55JnnhhdmleEMvAg4XzNYvcB4chQ5L4O+QgbOyBQ4ADQk4ADQk4ADQkIADQEMCDqyVZZ9nHi6KgEMj4nR3R+eZPzhIxrh5nnn/TqwjAYcmxOneVuE883BRBByaEKd7W4XzzMNFEXBoQpzubRXOMw8XRcDZaJ1+pyxO97Yq55mHiyDgbKxuv1MWp3tznnk2SY0xlr0O37S7uzuuXbu27NVgQ+zszKJ9u8uXZ9/OtYqm09nvvG/cmG15X70qTtBdVV0fY+ye+nECzqa6dGm25X27qtlXbAJchAcNuF3obCy/UwY6E3A2lt8pA50JOBvLAU9AZw8vewVgmSYTwQZ6sgUOAA0JOAA0JOAA99DpjH1sDr8DB7iLozP2HX2RzNEZ+xLHT7BctsAB7sK3wLGqBBzgLnwLHKtKwAHuwhn7WFUCDnAXztjHqhJwgLtwxj5WlaPQAe7BGftYRbbAAaAhAQeAhgQcABoScABoSMABoCEBB07kCzxgtfkzMuBFfIEHrD5b4MCL+AIPWH0CDryIL/CA1SfgwIv4Ag9YfQIOvIgv8IDVJ+DAi/gCD1h9jkIHTuQLPGC12QIHgIYEHAAaEnBgbTh7HJtEwIGlO4/wHp097uAgGePm2eNEnHUl4MBSnVd4nT2OTSPgwFKdV3idPY5NI+DAUp1XeJ09jk0j4MBSnVd4nT2OTSPgwFKdV3idPY5Nc6aAV9X3VtUXq+qFqtq97bYfr6qnqurJqnrX2VYTWFfnGd7JJHnmmeSFF2aX4s06O+upVJ9I8j1J/tXxwap6LMkHkrwpybcn+VRVvXGM8UdnfD1gDTltK5zembbAxxi/NcZ48oSb3pvkw2OMb4wxvpzkqSRvO8trAQA3Lep34K9J8pVjy8/OxwCAc3DPXehV9akkrzrhpitjjF+808NOGBt3eP69JHtJsu3vPQDgvtwz4GOMdz7A8z6b5HXHll+b5Lk7PP9+kv0k2d3dPTHyAMCtFrUL/WNJPlBVL6mqR5O8IclvLOi1AGDjnPXPyN5XVc8meXuSX6mqTyTJGOOLSf59kv+e5ONJfsQR6ABwfs70Z2RjjI8m+egdbruaxDmQAGABnIkNABoScABoSMABoCEBB4CGBBwAGhJwAGhIwAGgIQEHgIYEHAAaEnAAaEjAAaAhAQeAhgQcABoScABoSMABoCEBB1qZTpOdneTSpdnldLrsNYLleHjZKwBwv6bTZG8vOTycLR8czJaTZDJZ3nrBMtgCB9q4cuVmvI8cHs7Gz4OtezqxBQ60cePG6cZPw9Y93dgCB9rY3j7d+GkseusezpuAA21cvZpsbd06trU1Gz+rRW7dwyIIONDGZJLs7yeXLydVs8v9/fPZxb3IrXtYBAEHWplMkmeeSV54YXZ5Xr+fXuTWPSyCgANksVv3sAiOQgeYm0wEmz5sgQNAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4QFPTabKzk1y6NLucTpe9Rlykh5e9AgCc3nSa7O0lh4ez5YOD2XKSTCbLWy8uji1wgIauXLkZ7yOHh7NxNoOAAyzRg+4Gv3HjdOOsHwEHWJKj3eAHB8kYN3eD30/Et7dPN876EXCAJTnLbvCrV5OtrVvHtrZm42wGAQdYkrPsBp9Mkv395PLlpGp2ub/vALZN4ih0gCXZ3p7tNj9p/H5MJoK9yWyBAyyJ3eCchYADLInd4JyFXegAS2Q3OA/KFjgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADR0poBX1fdW1Rer6oWq2j02vlNVv19Vj89/fu7sqwoAHDnrFvgTSb4nyWdOuO1LY4w3z38+eMbXAWCFTafJzk5y6dLscjpd9hqtv4fP8uAxxm8lSVWdz9oA0M50muztJYeHs+WDg9lykkwmy1uvdbfI34E/WlWfq6pfr6p33OlOVbVXVdeq6trzzz+/wNUBYBGuXLkZ7yOHh7NxFueeW+BV9akkrzrhpitjjF+8w8N+O8n2GOPrVfXWJP+hqt40xvjd2+84xthPsp8ku7u74/5XHYBVcOPG6cY5H/cM+Bjjnad90jHGN5J8Y379elV9Kckbk1w79RoCsNK2t2e7zU8aZ3EWsgu9qh6pqofm11+f5A1Jnl7EawGwXFevJltbt45tbc3GWZyz/hnZ+6rq2SRvT/IrVfWJ+U1/KckXqurzST6S5INjjN8526oCsIomk2R/P7l8OamaXe7vO4Bt0WqM1fm18+7u7rh2zV52ADZHVV0fY+ze+563ciY2AGhIwAE2jJOurIczncgFgF6cdGV92AIH2CBOurI+BBxggzjpyvoQcIANcqeTqzjpSj8CDrBBnHRlfQg4wAZx0pX14Sh0gA0zmQj2OrAFDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4AGtvOk12dpJLl2aX0+my1+jsHl72CgDAIk2nyd5ecng4Wz44mC0nyWSyvPU6K1vgAKy1K1duxvvI4eFsvDMBB2Ct3bhxuvEuBByAtba9fbrxLgQcgLV29WqytXXr2NbWbLwzAQdgrU0myf5+cvlyUjW73N/vfQBb4ih0ADbAZNI/2LezBQ4ADQk4ADQk4ADQkIADQEMCDgANCTgANCTgANCQgANAQwIOAA0JOAA0JOAA0JCAA0BDAg4ADQk4ADQk4ADQUI0xlr0O31RVzyc5OOGmVyT52gWvziow781i3pvFvDfL3eZ9eYzxyGmfcKUCfidVdW2Msbvs9bho5r1ZzHuzmPdmWcS87UIHgIYEHAAa6hLw/WWvwJKY92Yx781i3pvl3Ofd4nfgAMCtumyBAwDHLCXgVfXuqnqyqp6qqh874fYPVtV/q6rHq+o/V9Vj8/G/VlXX57ddr6q/euwxn54/5+Pzn1de5JzuxxnmvVNVv39sbj937DFvnT/mqar6l1VVFzmn+3GGeU+Ozfnxqnqhqt48v639+33sfn+zqkZV7R4b+/H5456sqned9jmX5UHnvO6f7WP3u33ea/3ZPna/2+e91p/tqvrBqnr+2Bz+7rHbfqCq/tf85weOjZ/+/R5jXOhPkoeSfCnJ65N8S5LPJ3nstvt827Hr353k4/Prb0ny7fPrfybJV4/d79NJdi96Phc0750kT9zheX8jyduTVJL/mOQ7lz3X85r3bff5jiRPr9P7Pb/fS5N8Jslnj+aT5LH5/V+S5NH58zx0v8/ZdM5r/dm+y7zX+rN9p3nfdvvafbaT/GCSnznhsS9P8vT88mXz6y970Pd7GVvgb0vy1Bjj6THGHyT5cJL3Hr/DGON3jy1+a5IxH//cGOO5+fgXk/zxqnrJBazzeXjged9JVb06s/j9lzH7L+DfJvkb57vaZ3Ze8/6+JP9uYWt5/u4577mfSPKTSf7fsbH3JvnwGOMbY4wvJ3lq/nz3+5zL8sBzXvfP9txJ7/WJ1uWzPXevea/rZ/sk70ryyTHG74wx/k+STyZ594O+38sI+GuSfOXY8rPzsVtU1Y9U1Zcye+P/3gnP8/4knxtjfOPY2L+e7674xyu4u+ms8360qj5XVb9eVe849pzP3us5l+y83u+/lRd/yFu/31X1liSvG2P88n0+9r7+LZfoLHM+bu0+2/eY99p+tu/z/V67z/bc+6vqC1X1kap63T0e+0Dv9zICftKb8aItrjHGz44x/lSSf5jkH93yBFVvSvLPk/zwseHJGOM7krxj/vP957bG5+Ms8/7tJNtjjLck+ftJPlRV33a/z7lk5/F+//kkh2OMJ44Nt36/q+pSkp9K8g9O8dhVf7/PMuej+6zdZ/se817bz/Z9vt9r99me+6UkO2OMP5vkU0n+zT0e+0Dv9zIC/myS1x1bfm2S5+5w32S2e+KbuxKq6rVJPprkb48xvnQ0Psb46vzy95J8KLPdHKvkgec935X69fn165n9/uWN8+d87SmecxnO9H7PfSC3/T/0NXi/X5rZ73o/XVXPJPkLST42P8jnTo897b/lRTvLnNf5s33Hea/5Z/uu7/fcOn62M8b4+rE9SD+f5K33eOyDvd9LOADg4cx+cf9obh4A8Kbb7vOGY9e/K8m1+fU/Mb//+094zlfMr/+xJB9J8sGLntsC5/1Ikofm11+f5KtJXj5f/s3MPhhHBz68Z9lzPa95z5cvzf/jfv26vd+33f/TuXlg05ty60FsT2d24MypnrPZnNf6s32Xea/1Z/tO854vr+1nO8mrj11/X5LPzq+/PMmXMzuA7WXz6w/8fi/rH+A9Sf5nZv9v88p87J8m+e759Z/O7ECWx5P82tE/Tma7Vv/vfPzo55WZHfh0PckX5o/76aMPxSr9nGHe75+Pfz7Jf03yXceeczfJE/Pn/JnMT86zSj8POu/5bX/56D/+Y2Nr8X7fdt/b/8ftyvxxT+bY0agnPecq/TzonNf9s32Xea/1Z/tO854vr+1nO8k/O/a+/lqSP33ssT+U2YGpTyX5O2d5v52JDQAaciY2AGhIwAGgIQEHgIYEHAAaEnAAaEjAAaAhAQeAhgQcABr6/8tiLdpR2I+AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#check x and y in figure\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(8, 8)\n",
    "ax.scatter(x, y, c='blue', label='data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#split the last point\n",
    "x_last = x.values[-1]\n",
    "y_last = y.values[-1]\n",
    "x = x[:-1]\n",
    "y = y[:-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Case 1**\n",
    "\n",
    "使用MSE:\n",
    "\n",
    "*Mathod 1*\n",
    "\n",
    "$$ a = \\frac{n\\sum_{i}^{n}{y_i x_i}-\\sum_{i}^{n}{y_i}\\sum_{i}^{n}{x_i}}{n\\sum_{i}^{n}{x_i^2} - (\\sum_{i}^{n}{x_i})^2} $$\n",
    "\n",
    "$$ b = \\frac{\\sum_{i}^{n}{y_i} - a\\sum_{i}^{n}{x_i}}{n} $$\n",
    "\n",
    "*Method 2*\n",
    "\n",
    "$$ a = \\frac{\\bar{y_i x_i} - \\bar{y_i}\\bar{x_i}}{\\bar{x_i^2} - \\bar{x_i}^2} $$\n",
    "\n",
    "$$ b = \\bar{y_i} - a\\bar{x_i} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using method 1: \n",
      "a = -150.403\n",
      "b = 55.679\n"
     ]
    }
   ],
   "source": [
    "#method 1\n",
    "n = len(x)\n",
    "a = (n*np.sum(y*x) - np.sum(y)*np.sum(x)) / (n*np.sum(x**2) - np.sum(x)**2)\n",
    "b = (np.sum(y) - a*np.sum(x)) / n\n",
    "print('Using method 1: \\na = {:.3f}\\nb = {:.3f}'.format(a, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prediction: -9.104, label: -11.398\n"
     ]
    }
   ],
   "source": [
    "#predict\n",
    "def predict(x, a, b):\n",
    "    return a*x + b\n",
    "\n",
    "y_pred = predict(x_last, a, b)\n",
    "print('prediction: {:.3f}, label: {:.3f}'.format(y_pred, y_last))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2d340e4edd8>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#plot this lines\n",
    "x_seires = np.linspace(x.min(), x.max(), 100)\n",
    "y_l2_m1 = predict(x_seires, a, b)\n",
    "ax.plot(x_seires, y_l2_m1, color='red', label='MSE Method 1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using method 2: \n",
      "a = -150.403\n",
      "b = 55.679\n"
     ]
    }
   ],
   "source": [
    "#method 2\n",
    "a = (np.mean(y*x) - np.mean(y)*np.mean(x)) / (np.mean(x**2) - np.mean(x)**2)\n",
    "b = np.mean(y) - a*np.mean(x)\n",
    "print('Using method 2: \\na = {:.3f}\\nb = {:.3f}'.format(a, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prediction: -9.104, label: -11.398\n"
     ]
    }
   ],
   "source": [
    "#predict\n",
    "y_pred = predict(x_last, a, b)\n",
    "print('prediction: {:.3f}, label: {:.3f}'.format(y_pred, y_last))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2d340e591d0>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#plot this lines\n",
    "x_seires = np.linspace(x.min(), x.max(), 100)\n",
    "y_l2_m2 = predict(x_seires, a, b)\n",
    "ax.plot(x_seires, y_l2_m2, color='green', label='MSE Method 2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Case 2**\n",
    "\n",
    "使用MAE:\n",
    "\n",
    "$$ a = \\frac{\\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{y_i x_i / \\Lambda_i} - \\sum_{i}^{n}{x_i / \\Lambda_i}\\sum_{i}^{n}{y_i / \\Lambda_i}}{\\sum_{i}^{n}{1 / \\Lambda_i}\\sum_{i}^{n}{x_i^2 / \\Lambda_i} - (\\sum_{i}^{n}{x_i / \\Lambda_i})^2} $$\n",
    "\n",
    "$$ b = \\frac{\\sum_{i}^{n}{y_i / \\Lambda_i} - a\\sum_{i}^{n}{x_i / \\Lambda_i}}{\\sum_{i}^{n}{1 / \\Lambda_i}} $$\n",
    "\n",
    "$$ |y_i - ax_i - b| = \\Lambda_i $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At epoch 10, a = -145.949 and b = 53.761\n",
      "Early stopping at epoch 18, a = -145.408 and b = 53.527\n",
      "Finally, a = -145.410\n",
      "b = 53.528\n"
     ]
    }
   ],
   "source": [
    "#When using MAE, we need to loop several iterations\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "a_current = -10\n",
    "b_current = 0\n",
    "epoches = 1000\n",
    "eps = 1e-3\n",
    "logging_interval = 10\n",
    "for epoch in range(1, epoches+1):\n",
    "    factor = np.abs(y - a_current*x - b_current)\n",
    "    n = np.sum(1 / factor)\n",
    "    a_new = (n*np.sum(y*x/factor) - np.sum(y/factor)*np.sum(x/factor)) / (n*np.sum(x**2/factor) - np.sum(x/factor)**2)\n",
    "    b_new = (np.sum(y/factor) - a_new*np.sum(x/factor)) / n\n",
    "    if epoch%logging_interval == 0:#print the log\n",
    "        print('At epoch {}, a = {:.3f} and b = {:.3f}'.format(epoch, a_current, b_current))\n",
    "    if np.isclose(a_new, a_current, atol=eps) and np.isclose(b_new, b_current, atol=eps):#early stopping\n",
    "        print('Early stopping at epoch {}, a = {:.3f} and b = {:.3f}'.format(epoch, a_new, b_new))\n",
    "        break\n",
    "    else:\n",
    "        a_current = a_new\n",
    "        b_current = b_new\n",
    "print('Finally, a = {:.3f}\\nb = {:.3f}'.format(a_current, b_current))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prediction: -9.105, label: -11.398\n"
     ]
    }
   ],
   "source": [
    "#predict\n",
    "y_pred = predict(x_last, a_current, b_current)\n",
    "print('prediction: {:.3f}, label: {:.3f}'.format(y_pred, y_last))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2d3409149e8>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#plot this line\n",
    "y_l1 = predict(x_seires, a_current, b_current)\n",
    "ax.plot(x_seires, y_l1, color='black', label='MAE')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHSCAYAAAATyJnbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde3zO9f/H8cdnBzsYcy6SQyLnHEaEbQ5zyJyPGRKRSt8vv1BOc5gZOSZKohQTJWc52wx9oykKlUPmUJLUGNuF8fn9MS1ynO3a57q25/122+1qn8+u63pJ3z2+7+tzXZ+PYZomIiIi4lxcrB5ARERE0k4BFxERcUIKuIiIiBNSwEVERJyQAi4iIuKEFHAREREn5Gb1ADcqUKCAWaJECavHEBERyTS7d+/+wzTNgmm9n0MFvESJEsTGxlo9hoiISKYxDOPYg9xPL6GLiIg4IQVcRETECSngIiIiTsihjoGLiEjaXblyhZMnT2Kz2aweRe7C09OTokWL4u7uniGPp4CLiDi5kydPkitXLkqUKIFhGFaPI7dhmiZnz57l5MmTlCxZMkMeUy+hi4g4OZvNRv78+RVvB2YYBvnz58/QV0kUcBGRLEDxdnwZ/XekgIuISLoZhkG3bt1Sv09OTqZgwYIEBwcDcPr0aYKDg3nyyScpX748zzzzDABxcXF4eXlRpUqV1K+PP/74lscPDAykWLFimKaZuq1169b4+Pjcda74+Hjeeeed1O+jo6NTZ3oQd7r/2bNnqV+/Pj4+PvTr1++BHz8tdAxcRETSLWfOnOzbt4+kpCS8vLzYuHEjjzzySOr+0NBQgoKC+O9//wvAd999l7qvVKlS7Nmz557PkSdPHnbs2EHdunWJj4/n1KlT97zP3wF/+eWXH+BPdf88PT0JCwtj37597Nu3z67P9TetwEVEJEM0a9aMNWvWAPDJJ5/w7LPPpu47deoURYsWTf2+cuXKaX78zp07s2jRIgCWLl1K27Ztb9o/ceJEatSoQeXKlRk5ciQAb7zxBkeOHKFKlSoMGjQIgAsXLtC+fXvKli1LSEhI6qp+8+bNVK1alUqVKtGzZ08uXboEwLp16yhbtix169Zl6dKlt50tZ86c1K1bF09PzzT/uR6UVuAiIllJ//5wH6vZNKlSBaZNu+ePde7cmTFjxhAcHMx3331Hz5492bZtGwCvvPIKnTp1YsaMGTRq1Ijnn3+eIkWKAKQG9m9vv/029erVu+XxGzZsSO/evbl69SqLFi1i9uzZhIWFAbBhwwYOHTrErl27ME2Tli1bEhMTw/jx49m3b1/qCj86Oppvv/2W/fv3U6RIEerUqcOOHTvw8/OjR48ebN68mTJlytC9e3feffdd+vbtS+/evdmyZQuPP/44nTp1Sve/zoyigIuISIaoXLkycXFxfPLJJ6nHuP/WpEkTfv75Z9atW8fatWupWrVq6kvN9/sSuqurK3Xr1mXx4sUkJSVx48WvNmzYwIYNG6hatSqQsso+dOgQxYoVu+VxatasmfpqQJUqVYiLiyNXrlyULFmSMmXKAPDcc88xc+ZMAgMDKVmyJKVLlwaga9euzJ49O+3/cuxAARcRyUruY6VsTy1btmTgwIFER0dz9uzZm/bly5ePLl260KVLF4KDg4mJiaF69eppevzOnTvTpk0bRo0addN20zQZMmQIL7744k3b4+LibnkMDw+P1H92dXUlOTn5pjfH/ZujvsNfx8BFRCTD9OzZk9DQUCpVqnTT9i1btpCYmAhAQkICR44cue3q+F7q1avHkCFDbjq+Dikr/A8++IALFy4A8Msvv/D777+TK1cuEhIS7vm4ZcuWJS4ujsOHDwMwf/58AgICKFu2LEePHuXIkSNAyrF9R6EVuIiIZJiiRYumvtP8Rrt376Zfv364ublx7do1XnjhBWrUqEFcXNwtx8B79uzJf/7zn9s+vmEYDBw48JbtjRs35ocffqB27doA+Pj4sGDBAkqVKkWdOnWoWLEizZo1o3nz5rd9XE9PTz788EM6dOhAcnIyNWrUoG/fvnh4eDB79myaN29OgQIFqFu37h3fZV6iRAnOnz/P5cuXWb58ORs2bKB8+fL3/Hf2oIy7vWyQ2fz8/ExHuR54ZCQMGwbHj0OxYhAeDiEhVk8lInKrH374gXLlylk9htyH2/1dGYax2zRNv7Q+llbgtxEZCX36wPVXezh2LOV7UMRFRMQx6Bj4bQwb9k+8/5aYmLJdRETEESjgt3H8eNq2i4iIZDYF/Dbu9MbIB3jDpIiIiF0o4LcRHg7e3jdv8/ZO2S4iIuIIFPDbCAmB2bOheHEwjJTb2bP1BjYREXEcCvgdhIRAXBxcu5Zyq3iLiNxZdr+c6MaNG6levTqVKlWievXqbNmy5YGf437pY2QiIpJu2f1yogUKFGDVqlUUKVKEffv20aRJE3755Re7PqdW4CIikiGy8+VEq1atmnp1tQoVKmCz2VLvby9agYuIZCH91/Vnz28ZeznRKg9XYVpTXU70fi8n+vnnn1O1atWbLppiDwq4iIhkCF1OFPbv38/rr7/Ohg0b7vnnSS8FXEQkC7mflbI9ZefLiZ48eZI2bdrw8ccfU6pUqfu6T3roGLiIiGSY7Ho50fj4eJo3b05ERAR16tRJ85/rQSjgIiKSYe52OVE/Pz8qV65M7dq1Uy8nCv8cA//7a/r06Xd8/L8vJ1qgQIGbtjdu3JguXbpQu3ZtKlWqRPv27UlISCB//vyplxP9+01st3Pj5UQrVaqEi4sLffv2xdPTM/VyonXr1qV48eK3vf+MGTM4fPgwYWFhqX+O33///X7+lT0wXU5URMTJ6XKiziMjLyeqFbiIiIgTUsBFRESckAIuIiLihBRwERERJ6SAi4iIOCG7n8jFMIw4IAG4CiQ/yDvtRERE5GaZtQKvb5pmFcVbRCRrutflRP/WqlUrateufdO2UaNG8cgjj9z0WfD4+PhMmduZ6VSqIiKSbve6nCiknK3sm2++wcfHh6NHj1KyZMnUfQMGDGDgwIGZPbZTy4wVuAlsMAxjt2EYfTLh+URExAJ3u5wopFylq0WLFjddFlQeXGaswOuYpvmrYRiFgI2GYfxommbM3zuvR70P8EDnxbWnhLO/8r/NH9G44xCrRxERuS/9+/e/ryt7pUWVKlWYNi19lxOFlKiPHDmShx56iPbt2zNkyD+/W6dOncqCBQsAyJs3L1FRURn6Z8iK7L4CN03z1+u3vwPLgJr/2j/bNE0/0zT9ChYsaO9x0iR0UjDP7B/K4vdvPa+viIjc7G6XEz19+jSHDx+mbt26lClTBjc3t9TLiULKS+h79uxhz549ivd9susK3DCMnICLaZoJ1/+5MTDGns+ZkcIGfcE3Y8oSwnR4Hzr1fsvqkURE7up+Vsr2dKfLiS5evJi//vor9bj3+fPnWbRoEWPHjrVqVKdn7xX4Q8B2wzD2AruANaZprrPzc2YYn3wPsyb0R+qc86XLyelaiYuI3MOdLif6ySefsG7dOuLi4oiLi2P37t06Dp5Odg24aZo/m6b55PWvCqZphtvz+ezh74jXvR7xRbP/Y/VIIiIO63aXE42Li+P48ePUqlUrdVvJkiXJnTs3O3fuBFKOgd/4MbK4uLjMHNsp6XKi9+niX7/TfHQZtvmeI7JIPzq/+LbVI4mIALqcqDPR5UQtkDNvIdaMPEi9c76E/DqDT97rZ/VIIiKSjSngafB3xP3P5aHrrzNZOOsVq0cSEZFsSgG/i6+++orevXtz5cqV1G058xZi9cif8D+Xh26n3iHy3ZcsnFBERLIrBfwuvvrqK+bMmUOnTp24fPly6va/Ix5wLg/df5vFAkVcREQymQJ+F/379+ett95i2bJltGnThqSkpNR9OfMWYtWogwScy8NziriIiGQyBfwe/vOf//Dee++xdu1agoODuXjxYuq+nHkKsnr0IUVcREQynQJ+H/r06cO8efOIjo6madOmnD9/PnWft2+BlIifz0P307OY/05fCycVEbHeqFGjmDRp0h33L1++nAMHDmTiRFmTAn6funfvzieffMJXX31FUFAQf/31V+o+b98CrB51iPrn8vLc7+8xf+aLFk4qIuLYFPCMoYCnQceOHVmyZAl79uyhQYMGnDlzJnWft28BVo06SIP4vDx3ZjYfz9SVU0XEMUVGQokS4OKSchsZmf7HDA8P54knnqBRo0b89NNPALz//vvUqFGDJ598knbt2pGYmMiXX37JypUrGTRoEFWqVOHIkSO3/Tm5NwU8jVq1asXKlSv58ccfCQwM5NSpU6n7vH0LsHL0QRrG56PHmff5aGZvCycVEblVZCT06QPHjoFpptz26ZO+iP99XvNvv/2WpUuX8vXXXwPQtm1bvv76a/bu3Uu5cuWYO3cuTz/9NC1btmTixIns2bOHUqVK3fbn5N4U8AfQpEkT1q5dy7FjxwgICODEiROp+7x9C7Bi9E80jM/H82fmKOIi4lCGDYN/L3ATE1O2P6ht27bRpk0bvL29yZ07Ny1btgRg37591KtXj0qVKhEZGcn+/ftve//7/Tm5mQL+gAIDA9mwYQOnT5/G39+fo0ePpu7z9i3AyjGHUiM+b8YLFk4qIvKP48fTtv1+GYZxy7YePXowY8YMvv/+e0aOHInNZrvtfe/35+RmCng6PP3002zevJlz585Rr149Dh48mLrPK3c+Vo45RKP4fPT8Yy4fvt3TwklFRFIUK5a27ffD39+fZcuWkZSUREJCAqtWrQIgISGBwoULc+XKFSJveI0+V65cJCQkpH5/p5+Tu1PA08nPz4/o6GguX76Mv78/+/btS93nlTsfK8YcIig+P73OfsgH05+3cFIREQgPB2/vm7d5e6dsf1DVqlWjU6dOVKlShXbt2lGvXj0AwsLCeOqppwgKCqJs2bKpP9+5c2cmTpxI1apVOXLkyB1/Tu7BNE2H+apevbrprA4cOGAWLlzYzJ8/v/nNN9/ctC/x3FmzSf/8pjESc+5bPSya8PYWLDDN4sVN0zBSbhcssHoiEUmrAwcOpOnn9b9769zu7wqINR+gmVqBZ5By5coRExNDzpw5adCgQepF6iFlJb487DCNz+XnhT/nOcxK3B7vRhURxxcSAnFxcO1aym1IiNUTyYNQwDPQ448/TkxMDPnz56dRo0Zs27YtdZ+nT57UiPf6ax5zp/ewbtDr7PFuVBERyRwKeAYrXrw4W7dupWjRojRp0oRNmzal7vs74k3jC/DCXx8x563nLJzUfu9GFRER+1PA7eCRRx5h69atPP744wQHB7NmzZrUfZ4+eVgWdohm8QXpHf+xpRG3x7tRRcQaKYdSxZFl9N+RAm4nhQoVIioqiooVK9KmTRuWLl2aus/TJw9Lww6mRvz9ad0tmdEe70YVkczn6enJ2bNnFXEHZpomZ8+exdPTM8Me03Ckv3A/Pz8zNjbW6jEy1Llz52jWrBm7du3i448/pkuXLqn7bBfiaTeiDF/kOcNs32707v9xps8XGZlyzPv48ZSVd3i43tAi4myuXLnCyZMndQIUB+fp6UnRokVxd3e/abthGLtN0/RL6+Mp4JngwoULtGjRgq1btzJnzhx69vznpC43Rvw93xD69F9g4aQiIpLZHjTgegk9E/j4+LBmzRqCgoLo1asXM2fOTN3n6ZOHpWMP80x8QV48F8l7U7X8FRGRe1PAM4m3tzcrV66kZcuW9OvXj8mTJ6fu88iZOzXifc8vZNaULnd5JBEREQU8U3l4eLBkyRI6dOjAwIEDGTt27D/7rke8eXwhXkr4hHcnP2vhpCIi4ugU8Ezm7u7OwoUL6datGyNGjGDYsGGp7xz1yJmbz8ceonl8IV6+sEgRFxGRO3KzeoDsyM3NjXnz5uHl5cW4ceNISkpi8uTJGIaRGvEOw8vwMoswJ5u8/Noiq0cWEREHo4BbxMXFhVmzZuHp6cnUqVOx2WzMmDEDFxcXPHLm5rOxB+kwvAyvsBgmo4iLiMhNFHALGYbBtGnT8PLyYsKECdhsNt5//31cXV3xyJmbJeGH6DCsNK+wGHPSNV4Z+KnVI4uIiIPQMXCLGYZBREQEo0eP5sMPP6Rr165cuXIFgBzeufhs3GFaxT9Ev4ufMXNSR4unFRERR6GAOwDDMAgNDWXChAksWrSITp06cenSJQByePnw6bjDtIp/WBEXEZFUCrgDGTx4MNOnT2fZsmW0adOGpKQk4O+IH6L19YjPmNjB4klFRMRqCriDefXVV5k9ezbr1q0jODiYixcvAikRXzzuEK3jC/Nq4hLentje4klFRMRKCrgD6t27Nx999BHR0dE0adKE8+fPA39H/CBt4gvzn8TPmf5mO4snFRERqyjgDqpbt24sWrSInTt30qhRI/78808gJeKLrkf8v0lLeevNthZPKiIiVlDAHViHDh1YunQpe/fupUGDBpw5cwa4vhKPOEyb+ML0T1qmiIuIZEMKuINr0aIFq1at4uDBgwQGBnLq1CkA3D29b4r4tAltLJ5UREQykwLuBBo3bszatWs5duwY/v7+nDhxAvgn4m3PFWGAbbkiLiKSjSjgTiIgIICNGzfy+++/4+/vz9GjR4GUiC8ad4h25x5hgG05U8e3tnhSERHJDAq4E6lduzZbtmzh/Pnz1KtXj59++glIifgn4w7S/twj/N+lFUwZ38riSUVExN4UcCdTvXp1oqKiuHLlCgEBAezbtw9IifjC6xF/7dJKRVxEJItTwJ1Q5cqV2bp1K66urgQGBvLtt98C/0S8w7mivHZpJZMjWlg8qYiI2IsC7qTKli1LTEwMOXPmpH79+uzcuRNIiXjkuJ/ocK4oAy+vZlJEsMWTioiIPSjgTqxUqVJs27aNAgUK0KhRI2JiYoDrK/Hxh+h47lEGXV7DxHHNLZ5UREQymgLu5IoVK0ZMTAxFixaladOmbNq0CQC3HJ5Ejj9Ip3OPMvjKF4q4iEgWo4BnAUWKFGHr1q2ULl2a4OBg1qxZA6REfMENEX8z/BmLJxURkYyigGcRhQoVIioqiooVK9KmTRs+//xz4J+Idz5fjNeT1zIhvJnFk4qISEZQwLOQfPnysXnzZmrUqEGnTp1YuHAhkBLx+RE/0fl8Md5IXsf4sU0tnlRERNJLAc9ifH19Wb9+PfXq1aNr16588MEHwD8Rf/Z8cYZcXa+Ii4g4OQU8C/Lx8WHNmjUEBQXRq1cvZs6cCaRE/OOIH+lyPeIRY5tYPKmIiDwoBTyL8vb2ZuXKlbRs2ZJ+/foxadIkICXiH12P+NCrGxRxEREnpYBnYR4eHixZsoQOHTowaNAgxo4dC9wa8XFhjS2eVERE0koBz+Lc3d1ZuHAh3bp1Y8SIEQwbNgzTNFNeTh9/kJDzJRh2bSPhYUFWjyoiImngZvUAYn9ubm7MmzcPLy8vxo0bR2JiIlOmTMHVPQcfjf8J440nGJ57E4QFMWzERqvHFRGR+6CAZxMuLi7MmjULT09Ppk2bhs1mY+bMmbi652De+J/gjbIMz70JM6wRw0dssnpcERG5BwU8GzEMg2nTpuHl5cWECRNISkpi7ty51yP+I8YbZRmRezPmmIaMCN1s9bgiInIXCng2YxgGEREReHt7M3LkSGw2G/Pnz8fdPQcfjv8RY0g5QnNtwRzTgNDQLVaPKyIid6A3sWVDhmEQGhrKm2++yeLFi+nQoQOXLl3C1T0HH0T8QPeExxhpRjF6dH2rRxURkTtQwLOxQYMG8fbbb7NixQpat25NUlJSasSfSyjFKKIZNTrQ6jFFROQ2FPBsrl+/frz//vusX7+e4OBgLl68iKt7DuZGHKBHQilGs5VRowKtHlNERP5FARdeeOEFPv74Y6Kjo2nSpAnnzp3D1T0HcyIO0CPhcUYbiriIiKNRwAWArl27snjxYnbu3ElQUBB//vnn9Yjv5/kLpRVxEREHo4BLqvbt27N06VL27t1LgwYNOHPmTErExx+4IeIBYJpWj5pmkZFQogS4uKTcRkZaPZGISPoo4HKTFi1asGrVKg4ePEhgYCCnTp3CxdXthojHMHJ0IOa1a1aPet8iI6FPHzh2LOX/exw7lvK9Ii4izkwBl1s0btyYtWvXcuzYMfz9/Tlx4kRqxHteKM0YJ4v4sGGQmHjztsTElO0iIs5KAZfbCggIYOPGjZw5cwZ/f39+/vlnXFzdeH/8AXpdKEOYyzZCRwc4RcSPH0/bdhERZ6CAyx3Vrl2bzZs3c/78efz9/fnpp59wcXVj9vj9vHChDGNdtjtFxIsVS9t2ERFnoIDLXVWvXp3o6GiuXLmCv78/33//PS6ubrw3fj+9LzzBWJftjBjt79ARDw8Hb++bt3l7p2wXEXFWCrjcU6VKldi6dStubm7Ur1+fb775BhdXN2aN30fvi2UJd9nB8FH1HDbiISEwezYULw6GkXI7e3bKdhERZ2WYDvSRID8/PzM2NtbqMeQOjhw5QoMGDTh37hzr1q2jVq1aXLuaTN8hlXg/548Mvfo0Y0dtw3DR/y8UEblfhmHsNk3TL633029auW+lSpVi27ZtFChQgKCgIGJiYlJW4hHf0/tiWca5fsmwkXUddiUuIpKVKOCSJsWKFSMmJoaiRYvStGlTNm7cmBrxPhfLEeH2P4aOrKOIi4jYmQIuaVakSBG2bt1K6dKlCQ4OZvXq1bi4uvFuxHe8mFie8W5fMST0aUVcRMSO7B5wwzCaGobxk2EYhw3DeMPezyeZo1ChQkRFRVG5cmXatGnDkiVLcHF1451xe3kpsQIT3HfyRmhtRVxExE7sGnDDMFyBmUAzoDzwrGEY5e35nJJ58uXLx6ZNm6hZsyadOnUiMjISF1c3ZkZ8x0uJFXjTfZciLiJiJ/ZegdcEDpum+bNpmpeBRUArOz+nZCJfX1/Wr1+Pv78/3bp1Y+7cuRguLjdF/PURtRRxEZEMZu+APwKcuOH7k9e3pTIMo49hGLGGYcSeOXPGzuOIPfj4+PDFF1/QpEkTXnjhBWbOnJka8ZcTKzIxx9eKuIhIBrN3wI3bbLvpg+emac42TdPPNE2/ggUL2nkcsRcvLy+WL19Oq1at6NevH5MmTcJwcWFGxN7UiA8e8ZQiLiKSQewd8JPAozd8XxT41c7PKRbx8PDgs88+o2PHjgwaNIgxY8aAYTAjYi+vJFViUo5YBo2oqYiLiGQANzs//tdAacMwSgK/AJ2BLnZ+TrGQu7s7kZGReHp6MnLkSGw2G+Hh4bw9bg8MrcJkr90woiYTw3bpjG0iIulg14CbpplsGEY/YD3gCnxgmuZ+ez6nWM/NzY0PP/wQT09PIiIiSExMZOrUqbw9bg/G9Yibw2syaawiLiLyoOy9Asc0zS+AL+z9POJYXFxcmDVrFp6enrz11lvYbDbeeecdpl9fiU/x2o05vAaTx36tiIuIPAC7B1yyL8MwmDZtGt7e3owfPx6bzcbcuXOZPm4PxtCqTPX6BhRxEZEHooCLXRmGwbhx4/Dy8ko9Jj5//nzeGvctxrBqTPX8BnOYH1PCYxVxEZE0UMDF7gzDIDQ0FC8vLwYPHsylS5dYtGgR08K/wRhWnWme34IiLiKSJvptKZlm0KBBvP322yxfvpzWrVtju3SJqeG7+a+tCtM8v2XAsOr6iJmIyH1SwCVT9evXjzlz5rB+/XqaN2/OxcREpobvpr+tKm957lHERUTukwIuma5Xr17Mnz+fmJgYmjZtyvmEBKaExzLgUjXe8txD/2HVFHERkXtQwMUSISEhLF68mJ07dxIUFMRf8fFMHvs1/3epOtM99/LfoVUVcRGRu1DAxTLt2rVj6dKl7N27l/r163Pmjz+YNHYXr12uztte3/GfoVUUcRGRO1DAxVItWrRg9erVHDp0iMDAQE799hsTw1IiPsPre0VcROQOFHCxXFBQEGvXruXEiRMEBARw4uRJJobtYuBlP2Z4fc+ririIyC0UcHEIAQEBbNiwgd9//x1/f39+PnqUN8N2MvCyHzO9vqffkCcVcRGRGyjg4jBq167Nli1bSEhIwN/fn58OHuTNsJ0MulyDd7z38cqQyoq4iMh1Crg4lOrVqxMdHU1ycjIBAQHs27+fCWFfMfhKTd713s8rQypz7Wqy1WOKiFhOAReHU6lSJbZu3YqbmxuBgYF8u2cP48f8j8FXnuJd7/30G1pFEReRbE8BF4dUtmxZYmJiyJUrFw0aNOCrnTsZP+ZLXk+ulbISH/qkIi4i2ZoCLg6rVKlSxMTEULBgQRo3bkzMtm1EjN7BG8m1mOV9QBEXkWxNAReHVqxYMWJiYnj00Udp1qwZGzdtYtwNEX9Zx8RFJJtSwMXhFS5cmOjoaEqXLp1y4pc1axg3egdDkmvzXs4fFHERyZYUcHEKhQoVIioqisqVK9O2bVs+X7qU8NHbUyPed0glRVxEshUFXJxGvnz52LRpEzVr1qRTp05ELlxI+OjtDL36NO/n/JEX36ioiItItuFm9QAiaeHr68v69etp2bIl3bt3x2azMXbUNoxR/oT77IA3KvLe+H24uOo/bRHJ2vRbTpyOj48Pa9asoW3btvTu3RubzUbYqBiM0QGM9dmO+UYFZo/fr4iLSJam33DilLy8vFi+fDmdOnXi1VdfJSkpiTEjt2KMDiTMZxvmG+V5f/wBRVxEsiz9dhOn5eHhwWeffUa3bt0YPHgwNpuNUaFRMKY+YT7bQBEXkSxMv9nEqbm7uxMZGYmnpyehoaEkJSUxNiwKI6wBY3xiMN8ozxxFXESyIP1WE6fn6urKBx98gJeXFxERESQmJjJ1ajTG6PqM9tkKiriIZEH6jSZZgouLC++88w6enp5MmzYNm83GO+9sgTENFHERyZL020yyDMMwmDJlSupK3GazMXfuZghryGifrZivl2NOxH5c3XNYPaqISLop4JKlGIbBuHHj8Pb2ZsSIEdhsNubP34gxrjGjckXDkPLMiTigiIuI01PAJUsaPnw4np6eDBo0CJvNxuLF62B8U0blisYcUp65iriIODkFXLKsgQMH4uXlRb9+/WjdujVLl36B8WZzRuaKwhxSjg8iflDERcRpKeCSpb3yyit4enrSu3dvmjdvzsqVK2FKy9SIf6iIi4iTUt+3CskAACAASURBVMAly+vVqxeenp4899xzNG3alDVr1mC81ZbQXFtAERcRJ6WrkUm2EBISwuLFi9m5cyeNGjXilX6fMcZowPxcP9PjjbJcvXLZ6hFFRNJEAZdso127dixbtozvvvuO+vXr82LfTwhzaciC3Efp8cYTiriIOBUFXLKV4OBgVq9ezaFDhwgMDKRnr48Z69KIBbnjeE4RFxEnooBLthMUFMTatWs5ceIE/v7+dHtuLuNcGxOZO47ub5Qh+bLN6hFFRO5JAZdsKSAggI0bN/LHH3/g7+9Px2ffIcK1CQtzH+O5IWUVcRFxeAq4ZFu1atViy5YtXLhwAX9/f1q3n5Ya8e6KuIg4OAVcsrVq1aoRHR1NcnIyAQEBNG81kfFuTflEERcRB6eAS7ZXsWJFYmJicHd3JzAwkEZNx6ZGvNuQJxRxEXFICrgI8MQTTxATE0OuXLlo0KAB/g1CmeDWjEW5jyviIuKQFHCR6x577DFiYmIoVKgQQUFBPFX39dSId9W700XEwSjgIjcoVqwYW7dupVixYjRr1owqNfrzpvszLPY9oYiLiENRwEX+pUiRIkRHR1O6dGlatGhB2Up9UyMeooiLiINQwEVuo1ChQkRFRVG5cmXatm1LidI9mJijOZ/6nqDLG6UVcRGxnAIucgf58uVj06ZNPPXUU3Tu3JmHH+3MpBzBfOZ7ki5vlOaKLdHqEUUkG1PARe7C19eXdevWERgYSPfu3fEt2IrJHi35zPckIUOfUMRFxDIKuMg9+Pj4sHr1apo2bUrv3r3J4RPEFI9WKSvxoWUUcRGxhAIuch+8vLxYtmwZrVu35tVXXyXZtQ5TPVuzxPcXnlXERcQCCrjIffLw8ODTTz+lc+fODB48mPOXqjLFoxWfK+IiYgEFXCQN3N3dWbBgAT169GDkyJGcji+bGvHOQ/XGNhHJPAq4SBq5uroyd+5c+vbty4QJE4g7VZypnq1Z6vsrnRRxEckkCrjIA3BxceGdd96hf//+TJ8+nR+PPsRUj9Ys8/2VTkMe53LSBatHFJEsTgEXeUCGYTBlyhSGDBnCe++9x7c/5GZKjtYsy3OKTkPLKOIiYlcKuEg6GIbBuHHjCAsL4+OPP+arPTmY6tGG5XlO0XFYaUVcROxGARfJAMOHD2fixIl8+umnRH91jck52rDC9zc6DlXERcQ+FHCRDDJw4EBmzJjBihUr2LA1kclurVmRRxEXEftQwEUy0CuvvMKcOXPYsGEDqzbFM9mtDSvy/EYHRVxEMpgCLpLBevXqxYIFC9i2bRtL1vzGRJfWrFTERSSDKeAidtClSxcWL15MbGwsi1acYIKREvH2Qx/nUmKC1eOJSBaggIvYSbt27Vi6dCnff/89C5YcYYLZmlV5TtNhWGkuXTxv9Xgi4uQUcBE7Cg4OZvXq1Rw+fJgPF//IuGstWZXnNO2Hl1HERSRdFHAROwsKCmLdunWcPHmSuZH7GXulJavznKbdcK3EReTBKeAimcDf35+NGzfyxx9/MHv+HsIut2BNnt8VcRF5YAq4SCapVasWW7Zs4eLFi7zzYSyjkpqzJs/vtB3+OJcunLN6PBFxMgq4SCaqVq0a0dHRXL16lXc+/JrQi8/wRZ4ztB1eGtuFeKvHExEnooCLZLKKFSsSExODu7s7M+Z+xbDzzfgi7xnajiijiIvIfVPARSzwxBNPEBMTQ+7cuXl7zg5e/6spa/Ococ0IrcRF5P4o4CIWeeyxx4iJiaFQoULMmLONgWeCWJfnD0VcRO6LAi5ioUcffZSYmBiKFy/OjLnb6P9bI0VcRO6LAi5iscKFCxMdHU3ZsmV554MY+v3SgPW+f9B6xOOKuIjckQIu4gAKFizIli1bePLJJ5n1YQwvHgtkg+9ZRVxE7shuATcMY5RhGL8YhrHn+tcz9noukawgb968bNq0iVq1ajH74xh6HqmniIvIHdl7BT7VNM0q17++sPNziTi93Llzs27dOgIDA/kgcjvdfnyaDb5naTX8cZLO/2n1eCLiQPQSuoiDyZkzJ6tXr6Zp06Z8vPhLnt33FBvznKVVaGlFXERS2Tvg/QzD+M4wjA8Mw8h7ux8wDKOPYRixhmHEnjlzxs7jiDgHLy8vli1bRuvWrVn4+U7a76nBpjx/KuIiksowTfPB72wYm4CHb7NrGPAV8AdgAmFAYdM0e97t8fz8/MzY2NgHnkckq7ly5Qrdu3dn0aJFtGlejWXVv6HR+XysHH0Ir9z5rB5PRDKAYRi7TdP0S+v93NLzpKZpNrqfnzMM431gdXqeSyQ7cnd3Z8GCBXh6ejJv3jyaX36SNbX30jK0NCtG/4S3bwGrRxQRi9jzXeiFb/i2DbDPXs8lkpW5uroyd+5cXnrpJdZs3EtQTAU2+f5Jy5FlSDz3h9XjiYhF7HkM/E3DML43DOM7oD4wwI7PJZKlubi4MHPmTAYMGMDG6P0Ebi7LZt+/aKGIi2Rb6XoJ/W5M0+xmr8cWyY4Mw2Dy5Ml4eXkxbtw46lwpzZYmh2gxsgyrRh/Uy+ki2Yw+RibiRAzDIDw8nLCwMHbsOkTNNY8R5fMXwSNLayUuks0o4CJOaPjw4UyaNIld3/xMlZXFic4ZT/DI0lz863erRxORTKKAizip1157jZkzZ/Lt98eouOxRor3jaTH6CUVcJJtQwEWc2Msvv8zcuXPZ9+NJyn5ehGiveIIVcZFsQQEXcXI9e/ZkwYIFHDxymsc/fYitHoq4SHaggItkAV26dGHx4sXEnfiTEosKstU9nuajyyjiIlmYAi6SRbRr146lS5fy6+nzFP0kHzFu5xRxkSxMARfJQoKDg1m1ahV//JlE4ci8xLic45nRZbjw529WjyYiGUwBF8ligoKCWLduHecvXKHQAl+2medoPqasIn6DyEgoUQJcXFJuIyOtnkgk7RRwkSzI39+fjRs3YrsM+RfkYtvVczwz5glFnJRY9+kDx46Baabc9umjiIvzUcBFsqhatWoRFRWFaeQgz/yc7LhyXhEHhg2DxMSbtyUmpmwXcSYKuEgWVrVqVaKjo8nh6YPPfC922M7TLJtH/PjxtG0XcVQKuEgWV7FiRWJiYsiVOx/eCzz58mJKxBPO/mr1aJYoVixt20UclQIukg2UKVOGbdu2UaDgw3hE5uDL8+dpFlY2W0Y8PBy8vW/e5u2dsl3EmSjgItlEyZIliYmJoWjR4rgvdOd/fyVky4iHhMDs2VC8OBhGyu3s2SnbRZyJYZqm1TOk8vPzM2NjY60eQyRLO3XqFI0aNeLgwZ+41vEqFTxzcTn6R0aMKaKIiVjAMIzdpmn6pfV+WoGLWMDKzyEXLlyYfv2iuXq1EuYnruy/mECOwLIMeOWkPkol4kQUcJFM5gifQ54woSCmuQXTrI652IX95xJ4qH05xoaezLwhRCRdFHCRTOYIn0NO+chUXmAjJk9z7XODA39cwK1+Oc6fUcRFnIECLpLJHOFzyP98ZCo3sA6oz7UVcOC3CzQZp4iLOAMFXCSTOcLnkG/+KFVOYDUuLs9wbQ3s+iEl4ud+15lNRByZAi6SyRzhc8i3fpTKiw8+WEabNm24th52fX+BphEVFHERB6aAi2QyR/kcckgIxMXBtWspt889l4PFixfTuXNnrm2Gnd9coLFW4iIOS58DF5GbXL16lRdeeIF58+Zh1IYaNb3ZMPQHfAvpXKMi9qDPgYtIhnB1dWXu3Lm89NJLmP+D2C8TCRpXjvjTx6weTURuoICLyC1cXFyYOXMmAwYM4NrXEBuTSONx5RVxEQeigIvIbRmGweTJkxk6dCjmtxAblUiQIi7iMBRwEbkjwzAIDw8nLCwM83vYvTGRRuGKuIgjUMBF5J6GDx/OpEmTMH+Ab9Yl0nBsOf46ddTqsW7LyvPMi2QmBVzEiVgZp9dee42ZM2diHoJvv0iiYXgFh4u4I5xnXiSzKOAiTsIR4vTyyy8zd+5cOAp7VifRYKxjRdwRzjMvklkUcBEn4Shx6tmzJ5GRCzFOGOxdmUSDMMeJuCOcZ14ksyjgIk7CkeL07LPP8tlnS3A95cLe5UkEji7vEBF3hPPMi2QWBVyyNWd6w5Ojxalt27YsX7EStz9c+X6ZjYCR5fnz1yPWDHOdI5xnXiSzKOCSbTnCMeW0cMQ4NW/enLVr15PjnDvfL7URMLKipRF3lPPMi2QGnQtdsq0SJVKi/W/Fi6dc3MMRRUamHPM+fjxl5R0e7hhx2rZtG02aBJHkcYkKbT2ICdtPviKlrB5LxCk86LnQFXDJtlxcUlbe/2YYKVfokrTZuXMnjRoEcsHNRvm2HmwLV8RF7ocuZiKSRo52TNnZPfXUU8Rs/5Lc17w5sOQSdYdUsPyYuEhWpoBLtuWIx5SdXdWqVfnyq13kdfXhhyWXeHpwBc6ePGT1WCJZkgIu2Zbe8GQfFSpU4KtduyngmZufll7i6cGVFHERO9AxcBGxi6NHj1L7qWqcPh/P461y8L9J31Hg0SesHkvE4egYuIg4lJIlS/L17u94JF9+Di+7zFMDKvHHiZ+sHksky1DARcRuHn30UWK/2Ufxhwrx84or1PyPIi6SURRwEbGrhx9+mNhv91OqaGGOrrqC3yvOF3FnOmOfZB8KuIjYXYECBfj6m/2ULfkox9ZcofpLlThz/Aerx7ovznbGPsk+FHARyRR58+Zl5+59VCxTguNfXKFanyedIuKOchU4kX9TwEUk0+TOnZuvYvdRteLjnFx/hao9Kzt8xB3pKnAiN1LARSRT5cyZkx07v6NmtbL8sjmZJ59z7JfTdcY+cVQKuIhkOi8vL7b9by91a1bkVPRVKoVU4ve4/VaPdVs6Y584KgVcRCyRI0cOtmz/hvp1qnB6+1UqPvskp4/us3qsW+iMfeKodCY2EbHU1atXCW5cm3Vbvia/nyv7P93DQyUrWj2WSKbRmdhExCm5urqyZuNXtGzyNGdjr1K+bRV++/l7q8cScXgKuIhYzsXFheVrt9OhRQB/7rlKuVZV+PXwXqvHEnFoCriIOATDMFi8Ioqu7YOI33eN8i2q8cuhPVaPJeKwFHARcRiGYTD/sw30CnmGcz9eo1zzahz/cbfVY4k4JAVcRG7LyvN/z1mwhpefb0XCIZMKzWsQ98PXmffkIk5CAReRWzjC+b9nfrCcAX07cOFnkwpNn+Ln/Tsz78lFnIACLiK3cJTzf09591Nef7ULicdNKjatzeHv/pe5A4g4MAVcRG7hSOf/Hj89khGvPUfSLyaVm9bhx2+3Zf4QIg5IAReRWzja+b/HTJpH2Ou9SDptUqVpAPtjo60ZRMSBKOAicgtHPP/38Ig5TBjel0tnTao1a8B3OzdbN4yIA1DAReQWjnr+78Gj32Xq6Fe5fM6kxjON+GbHBmsHErGQzoUuIk5n5oT/o1/oVNy9YPvy1dQMbG71SCIPTOdCF5Fs45XXpzA7YjBXLkGdlsHs2LTC6pFEMp0CLiJOqff/TeDDN4eQnAwBrVuzde3nVo8kkqkUcBFxWj1eHcf8yaFcAxq2a0/o4IWWnT1OJLO5WT2AiEh6dH1pNG5urnR5bSTh00PI45qMaXZPPXscWP/mOxF70ApcRCyX3vOud+4dit8jYzG94K/k58jrNRew5uxxIplFARcRS2XUeddjfxpGac/xkBPik18gr/e7gDVnjxPJDAq4iFgqo867XqwYHPztdUp7vonhC/GXXiav93TLzh4nYm8KuIhYKqPOu/732eMOnhrE4+5TMPLDX7b/8kzDiekfUsQBKeAiYqmMOu/6jWePO/TbAJ4qMA2Xh+DdeYOZOWVs+gcVcTAKuIhYKiPPux4SAnFxcO0afLn/v2yc8TauRaDfwBFMiQjNkHlFHEW6Am4YRgfDMPYbhnHNMAy/f+0bYhjGYcMwfjIMo0n6xhSRrMqe511v0LYfm2bMxLU4vDY0jPGjXk//g4o4iPR+Dnwf0BZ478aNhmGUBzoDFYAiwCbDMMqYpnk1nc8nIllQSIj9Pqsd2OplolzdaND/RYaMfhPbpUuMiphmnycTyUTpWoGbpvmDaZo/3WZXK2CRaZqXTNM8ChwGaqbnuUREHlS94D5snT4HtzIwevxbvDHgJatHEkk3ex0DfwQ4ccP3J69vu4VhGH0Mw4g1DCP2zJkzdhpHRLK7p5/pRcxbc3EvDxOmzeK/L/XEka7GKJJW9wy4YRibDMPYd5uvVne722223fZ/KaZpzjZN0880Tb+CBQve79wiImlWu2lPYiZ/gHslmD7rQ17q1VURF6d1z2Pgpmk2eoDHPQk8esP3RYFfH+BxREQyVK2mz7PdxaDeG8/z3ocLsV26xAfzP8XFRR/KEedir/9iVwKdDcPwMAyjJFAa2GWn5xIRSZOajXuwfcI8PGrARws/J6Rja65e1Xtsxbmk92NkbQzDOAnUBtYYhrEewDTN/cCnwAFgHfCK3oEuIo6kRtBzbAufh8fTsOjzVbRr1Yzk5GSrxxK5b4YjHf/x8/MzY2NjrR5DRLKR3ZsXUHdMN2wx0CwokOWr15MjRw6rx5JsxDCM3aZp+t37J2+mgz4ikq1Vb9iV7aHz8WwIazdG80yTBthsNqvHErknBVxEsr3qDbuyfch8vJvA5ugdBDXwJ/Hfl0gTcTAKuIgIKRHfNngB3s/A9q++pkFAHRISEqweS+SOFHARkeuqNQhh22sLyNkCdu7eQ0DdWsTHx1s9lshtKeAiIjeo1iCEbQMW4tMavt13gLpP1+Ts2bNWjyVyCwVcRORfqgY+y7b/LsKnLew/eIina9Xg9OnTVo8lchMFXETkNqoEdGJbv0Xkag+H4o7ydK0a/PLLL1aPJZJKARcRuYMqAZ2IeWkRPp3g6K8nqF2rBseOHbN6LBFAARcRuasqAZ3Y1mcxuTrByTOnqF2rJocPH7Z6LBEFXETkXp7070jMC4vJ/SycPvc7dWrX4sCBA1aPJdmcAi4ich+e9O/I1ucXk+tZ+CPpLPXq1mHv3r1WjyXZmAIuInKfnvTvSEyPz8jdBc4lxxPgXw9dv0GsooCLiKRB5Xrt2drtM3KHwAUjgfqBAezYscPqsSQbUsBFRNKocr32bA1ZQu4uYMuRSOOgRkRFRVk9lmQzCriIOJXISChRAlxcUm4jI62Zo1LddmztkhLxyzltNGvWlHXr1lkzjGRLCriIOI3ISOjTB44dA9NMue3Tx9qIR3deQu5n4ZrvZVq2aMHy5cutGUayHQVcRJzGsGHw76t8JiambM8ID7K6r1S3HVs7LSV3Z6BgMu3bt2Px4sUZM5DIXSjgIuI0jh9P2/a0SM/qvmKdNkR3XErujuBa+BpdunTho48+Sv9QInehgIuI0yhWLG3b0yK9q/uKddqwteNyfNuD26PX6NGjB++99176BxO5AwVcRJxGeDh4e9+8zds7ZXt6ZcTqvsLTrYjqsBzfduDxGPTt25dp06alfziR21DARcRphITA7NlQvDgYRsrt7Nkp29Mro1b3FZ5uRXSHFeRuAx5lYMCAAURERKR/QJF/UcBFxKmEhEBcHFy7lnKbEfGGjF3dl6/Vkuj2K8jdEjzLwdChQwkNDcU0zYwZVgQFXEQEyPjVfflaLYlut4LczQ28KkFYWBivv/66Ii4Zxs3qAUREHEVISMat6OF6xI2VBBotwc1k4sSJJCUl8dZbb+HiovWTpI8CLiJiR+WeCiaalQTSEsPNZMaMGdhsNmbNmoWrq6vV44kTU8BFROzs3xGfM2cOSUlJzJs3Dzc3/RqWB6P/ckREMkG5p4LZaqymPi3A/RqRkZHYbDYWLlxIjhw5rB5PnJAOwoiIZJKyNZ8hqv0qcj3lgk8D+Pzzz2nXrh02m83q0cQJKeAiIpmobM1niO6wmlzVXcgVBKtXr6Zly5Yk/vs0cCL3oICLiGSyJ2o0I6rDanyquJCrGWzevJlmzZqRkJBg9WjiRBRwERELpEa8ogs5m19jx47tNG7cmPj4eKtHEyehgIuIWOSJGs2I7vgFucq6kLPlNWJjY2nYsCFnz561ejRxAgq4iIiFyvg1IarjF/g87oJ362T279tHYGAgp0+ftno0cXAKuIiIxcr4NSG60zp8Srrg0fYyRw4fIiAggF9++cXq0cSBKeAiIg6gdPWglIg/6kKOdpc4eeI4/v7+xMXFWT2aOCgFXETEQaRGvLALbu2T+OPMafz9/Tl06JDVo4kDUsBFRBxI6epBRHfZgE8hF4z2F7lwPh5/f38OHDhg9WjiYBRwEREH83jVhkR32UCu/K5ca59A8mUbAQEB7Nmzx+rRxIEo4CIiDujxqg2JDtlArryuJLeLx83FpH79+uzatcvq0cRBKOAiIg6qVJUGRIdswMfXFVvbs+T0ykGjRo3Yvn271aOJA1DARUQc2N8Rz+XjysXWv5Mvjw9NmjRhy5YtVo8mFlPARUQcXKkqDYgK2UCunK6ca3mKwoXy0bx5c9auXWv1aGIhBVxExAmUqtKA6K6byOXlytkWJynx6MO0aNGKQoWW4+ICJUpAZKTVU0pmUsBFRJzEY08GpkTcw5XfnonDw7U4Z860xzQXcewY9OmjiGcnCriIiBP5O+I5DVdyvHQYL/fSQAjwEYmJMGyY1RNKZlHARUQsFBmZ8vJ3Wl4Gf+zJQIx5W1Ii/uqPeOeoAPQA3uP4cbuOKw5EARcRsUhkZMrL3seOgWmSppfBXX39Yd4WfK654v7q9+T0qAr0JU+eaXafWxyDAi4iYpFhwyAx8eZt9/syeHg4/HXZH+ZF43PVFbdXv8XXuyZ//TWAiIgI+wwsDkUBFxGxyJ1e7r6fl8FDQmD2bHDzrQvzosmV7Ar9dtGs4dMMHTqU0NBQTNPM2IHFoSjgIiIWKVYsbdv/LSQE4uLg5F912dYrmrxX3fiy+pe0aupPWFgYgwcPVsSzMAVcRMQi4eHg7X3zNm/vlO1pVaJiXaKfT4l49JMxdGjZkEmTJvHqq69y7dq1jBlYHIoCLiJikb9fBi9eHAwj5Xb27JTtD6J4hTopEU92Y0O5zXTt0IyZM2fSp08frl69mrHDi+UMR3p5xc/Pz4yNjbV6DBERp3Zs/w7qfxjIWbdkOv7SijkLVtClSxc++ugj3NzcrB5P/sUwjN2mafql9X5agYuIZDHFK9Qh6vlo8ie78VnRFbz8fHsWLlxI586duXz5stXjSQZRwEVEsqDiFeoQ3TOGfFfciHxoCQP6Psvnn39Ou3btsNlsVo8nGUABFxHJooqVr010zxjyX3Hjgzyf8Pp/urN69WpatmxJ4r8/gC5ORwEXEcnCipWvTdT1iL+b82NCB/Zi8+bNNGvWjISEBKvHk3RQwEVEsrhi5WsT3Ws7BS67Mc19LqNff4EdO3bQuHFj4uPjrR5PHpACLiKSDTxa7imie22n4CU3JjKbiOF92b17Nw0aNOCPP/6wejx5AAq4iEg28Wi5p4ju8yUFL7kx9vJMJox4iR9++IH69evz22+/WT2epJECLiKSjRR9ogbRfb6k0CV3RiZNZ+LIl/n5558JCAjg5MmTVo8naaCAi4hkM0WfqEH0i1/y0CV3hp6fwpTRr3Dq1Cn8/f05evSo1ePJfVLARUSyoUfK+KVGfNCfE5k65hXi4+Px9/fn0KFDVo8n90EBFxHJpv6O+MOX3Bnw+3imjnkFm82Gv78/Bw4csHo8uQcFXEQkG3ukjB9R1yPe79RYpoW9DEBAQAB79uyxeDq5GwVcRCSbe6SMH9F9v6KwzZ2+J8YwbcxLeHl5Ub9+fXbt2mX1eHIHCriIiFCkdLXUiL8QN5Ipo/qQL18+GjVqxPbt260eT25DARcREeCfiBexufP8kRFMCu1FkSJFaNKkCVu2bLF6PPkXBVxERFIVKV2NqOsR735wGG8O60GpUqVo3rw5a9eutXo8uYECLiIiNylSuhrRL++iiM2dLj8OIeKNrpQvX55WrVqxbNkyq8eT6xRwERG5ReFSVYh+eRdFbTnotP91wgZ2onr16nTo0IFFixZZPZ6ggIuIyB0ULlWFqJd38mhSDjrue52R/21LnTp1CAkJYd68eVaPl+0p4CIickeFS1Uhqt/XPJqUg/bfD2ZYv5Y0bNiQ559/nlmzZlk9XraWroAbhtHBMIz9hmFcMwzD74btJQzDSDIMY8/1L/0ti4g4qYcfq5wa8bZ7BjKoT1OCg4N56aWXmDZtmtXjZVtu6bz/PqAt8N5t9h0xTbNKOh9fREQcwN8RbzCjBq2/fY2VPabg4eHBgAEDSExMpHjxoQwbBsePQ7FiEB4OISFWT521pSvgpmn+AGAYRsZMIyIiDuvhxyoT9epu6r9dnZbf/B8ru0zG09OTYcOG4eaWRHLyGMDg2DHo0yflPoq4/djzGHhJwzC+NQxjq2EY9ez4PCIikkkeKlmRqFd3UyLJg5a7X+P5NtXw8XmB5OSxwGDABCAxEYYNs3TULO+eATcMY5NhGPtu89XqLnc7BRQzTbMq8H/AQsMwct/h8fsYhhFrGEbsmTNnHuxPISIimeahkhXZ8mpsSsS/eY3H8pcD+gGTgFeBa0DKy+liP/cMuGmajUzTrHibrxV3uc8l0zTPXv/n3cARoMwdfna2aZp+pmn6FSxY8EH/HCIikolujPihrq9RufhjwCBgJtAHuEqxYtbOmNXZ5SV0wzAKGobhev2fHwNKAz/b47lERMQaD5WsSNR/dlP8ggeHuv4flYs/BIQCc3F17c6YMclWj5ilpfdjZG0MwzgJ1AbWGIax/vouf+A7wzD2AkuAvqZp/pm+UUVExNEUKlGBrf13UzzBg0NdB1KlZG7y5Ing6tWFrFzZmcuXL1s9YpZlmKZp9Qyp/Pz8zNjYWKvHEBGRNPo9bj8Np1fniPclVlV9k30nc9C/f3+aN2/OkiVL8PT0/P/27j64qvrO4/j7RwokKQQJT5KGJJhSC0FnnaYP7opiVp6UCAhotPWRIAAADL9JREFUlamtnR22zs6sohWJlK7Io01ZKMvOBBja4G4gjsEs2xISSXjIuMJYaCLyEOTBgKHYFReJBFKZ9rd/3JP0Jk1ISG7uuefcz2vmzj333HPO/X7zy49P7r1nDm6XGLGMMYestZkdb9mSrsQmIiLdNjQtg93PVpF+tS9Tq+YzNvkL8vLy2LFjB9nZ2TQ0NLhdou8owEVEokxBAaSlQa9egfuCgtAcd0jKaHY/W8VXnRD/6uB68vPz2b17N1OmTKG+vj40LySAAlxEJKoUFAQusnL2LFhL80VXQh3io67GMrVqPl/58h/YsmUL77zzDhMnTuTSpUuheSFRgIuIRJOFCwMXWQkW6ouuDEkZTcWzv2PU1Viyq19kUMxZtm3bRlVVFVlZWVy8eDF0LxbFFOAiIlGkvYurhPqiK0NSRrN7XnVziMd/cZzt27dTU1PD+PHj+fjjj0P7glFIAS4iEkXau7hKT1x0ZfCI29k9r5qvNcTy8Hs59KqvoqSkhNraWu677z7q6upC/6JRRAEuIhJFli2D+PiW6+LjA+t7wuARt1PxXCDEpx1+ieufHKCsrIwLFy5w7733Ultb2zMvHAUU4CIiUWTOHNiwAVJTwZjA/YYNPfu/hjWF+O0NsTz8/ks0nK+koqKCzz77jHHjxnHy5Mmee3Ef04VcREQkLD6tO8nfr7qTmn6NbB+7lFu/PpUJEyYQExNDeXk5GRkZbpfoCl3IRUREItqg5FFUPH+Y0VfimHbkJ1w4/mv27t2LMYbx48dTXV3tdomeogAXEZGwGZQ8ivLn32PMlTimH13EucNvsm/fPuLi4rj//vt599133S7RMxTgIiISVoOSR1H+wvvNIX760OtUVlaSmJjIAw88wNtvv+12iZ6gABcRkbBLTEpvDvFpRxdxfP9rVFZWkpSUxKRJk6ioqHC7xIinABcREVc0hfjYK3FMP/YvHK7cxL59+0hPT+ehhx6ipKTE7RIjmgJcRERck5iUzq6gED9UkceePXvIyMhg+vTpFBcXu11ixFKAi4iIq4JDfMbxlzlQ+m9UVFSQmZnJ7NmzKSwsdLvEiKQAFxER1yUmpVM+/yhjr8TxSM1i/mfHGsrKyrjnnnt44oknyM/Pd7vEiKMAFxGRiDBw+EjK5x/ljivxPFKzmH3bcykpKWHChAk89dRT5OXluV1iRFGAi4hIxBg4fCS75h/hzivxzDyxhN1vrmT79u1kZ2fz9NNPs2bNGrdLjBgKcBERiSgDh4/kLSfEH/lgCeVFyykqKmLWrFnMmzeP5cuXu11iRPiS2wWIiIi0NnD4SHYtOMaElWN45IMlbHvdsnXrVmJjY1m4cCHXrl3jlVdewRjjdqmuUYCLiEhEumVYKrsWHGPiygxmnlxK0dY/k5+fT2xsLEuXLuXatWvk5uZGbYgrwEVEJGLdMiyVtxYcZdLKDGadXM4bBZb169cTGxvLqlWruHr1KuvWraNXr+j7RlgBLiIiEe2WYam8lXOMiSvGMPvUCt74T1i7di1xcXHk5ubS2NjIxo0biYmJcbvUsFKAi4hIxBswNIW3co4xacUYZp1aQdF/wKuvvkp8fDyLFy+msbGRzZs307t3b7dLDRsFuIiIeMKAoSmUBYX4G69ZXn55BbGxseTk5NDY2EhhYSF9+vRxu9SwiL4vDURExLOaQjzz837MPr2S/9qcw4IFC1izZg3FxcXMmDGDxsZGt8sMCwW4iIh4yoChKZTmHG0O8eL8F3nmmWdYv349O3fuZOrUqTQ0NLhdZo9TgIuIiOcMGJpC2UvHyfy8H4+e+RnF+S8yd+5c8vPz2bNnD1OmTKG+vt7tMnuUAlxERDwpYUhyixB/81fzefLJJ9m6dSv79+9n4sSJXLp0CYCCAkhLg169AvcFBa6WHhIKcBER8aymEP9mfT8e+zCXbb98gUcffZSioiKqqqrIysoiL+8ic+fC2bNgbeB+7lzvh7ix1rpdQ7PMzEx78OBBt8sQERGPqf+kjsnLRvPbhCsUpv2YmT/MpbS0lBkzZvCnP6Vz/Xo5cGuLfVJTobbWlXJbMMYcstZm3ux+egcuIiKelzAkmdKFx/lWfX8eO/tztv3yBSZPnsyOHTu4fr0WuBf4qMU+5865UWnoKMBFRMQXEoYkU7qohm9fDoR40abnycrKYtiwMuBjAiH+YfP2KSluVRoaCnAREfGN/oOSKF1Uw3cu9+e75/6VNzY9x6pVf0ffvhXAZeDnAMTHw7JlrpbabQpwERHxlf6Dkti5qIa7Lyfw+LnV9Gl8jk2bvklS0gFgNampsGEDzJnjdqXdo0upioiI7/QflETJouM8uGQ0j59bzZYRlvPnV7tdVkjpHbiIiPhSU4jffTmBJz5aw+sbn3G7pJBSgIuIiG/1H5TEzp+e4G8vJzCnbq2vQlwBLiIivtYv8VZKfnqCey4PoL7hktvlhIy+AxcREd/rl3grFbn/S0xv//xXo3oHLiIiUcFP4Q0KcBEREU9SgIuIiHiQAlxERMSDFOAiIiIepAAXERHxIAW4iIiIBynARUREPEgBLiIi4kEKcBEREQ9SgIuIiHiQAlxERMSDFOAiIiIepAAXERHxIAW4iIiIBynARUREPEgBLiIi4kEKcBEREQ8y1lq3a2hmjPkEONvGU4OBi2EuJxKo7+iivqOL+o4uN+o71Vo75GYPGFEB3h5jzEFrbabbdYSb+o4u6ju6qO/o0hN96yN0ERERD1KAi4iIeJBXAnyD2wW4RH1HF/UdXdR3dAl53574DlxERERa8so7cBEREQniSoAbYyYbY04YY04ZYxa08fyPjDHvG2OqjTFvG2PGOOsnGGMOOc8dMsZkBe2z1zlmtXMbGs6eOqMbfacZY64F9ZYXtM83nH1OGWPWGmNMOHvqjG70PSeo52pjzJ+NMX/jPOf58Q7abpYxxhpjMoPW5Tj7nTDGTLrZY7qlqz37fW4Hbde6b1/P7aDtWvft67ltjPmBMeaToB7+Iei57xtjTjq37wetv/nxttaG9QbEAKeB24A+wHvAmFbbJAQtPwyUOst3AUnO8ljgfNB2e4HMcPcTpr7TgCPtHPdd4G7AADuBKW73Gqq+W21zB3DGT+PtbNcfqAQONPUDjHG27wuMdI4T09ljerRnX8/tG/Tt67ndXt+tnvfd3AZ+AKxrY99E4IxzP9BZHtjV8XbjHfi3gFPW2jPW2i+AQmBa8AbW2vqgh18GrLO+ylr7e2f9USDWGNM3DDWHQpf7bo8xZjiB8NtvA78BrwHTQ1t2t4Wq78eBrT1WZeh12LdjCfAzoDFo3TSg0Fr7R2vth8Ap53idPaZbutyz3+e2o62xbpNf5rajo779OrfbMgnYZa39P2vtJWAXMLmr4+1GgH8F+CjocZ2zrgVjzD8ZY04TGPh/buM4M4Eqa+0fg9b9yvm4YlEEftzU3b5HGmOqjDH7jDHjgo5Z19ExXRaq8X6Mv57knh5vY8xdwAhr7W86uW+nfpYu6k7PwXw3tzvo27dzu5Pj7bu57ZhpjDlsjCkyxozoYN8ujbcbAd7WYPzVOy5r7b9ba9OBF4GftDiAMRnAq8A/Bq2eY629Axjn3L4XsopDozt9XwBSrLV3Ac8BW4wxCZ09pstCMd7fBq5aa48Erfb0eBtjegGrgedvYt9IH+/u9Ny0je/mdgd9+3Zud3K8fTe3Hb8G0qy1dwLlwOYO9u3SeLsR4HXAiKDHycDv29kWAh9PNH+UYIxJBoqBJ621p5vWW2vPO/efA1sIfMwRSbrct/NR6qfO8iEC3798zTlm8k0c0w3dGm/Hd2n1F7oPxrs/ge969xpjaoHvAP/tnOTT3r43+7MMt+707Oe53W7fPp/bNxxvhx/nNtbaT4M+QdoIfKODfbs23i6cAPAlAl/cj+QvJwBktNpmVNByNnDQWb7F2X5mG8cc7Cz3BoqAH4W7tx7sewgQ4yzfBpwHEp3HvyUwMZpOfHjQ7V5D1bfzuJfzy32b38a71fZ7+cuJTRm0PIntDIETZ27qmB7r2ddz+wZ9+3put9e389i3cxsYHrQ8AzjgLCcCHxI4gW2gs9zl8XbrB/Ag8AGBvzYXOuteAR52ln9B4ESWamBP0w+HwEerDc76pttQAic+HQIOO/v9omlSRNKtG33PdNa/B/wOyA46ZiZwxDnmOpyL80TSrat9O8+Nb/rlD1rni/FutW3rf9wWOvudIOhs1LaOGUm3rvbs97l9g759Pbfb69t57Nu5DawIGtc9wNeD9v0hgRNTTwFPdWe8dSU2ERERD9KV2ERERDxIAS4iIuJBCnAREREPUoCLiIh4kAJcRETEgxTgIiIiHqQAFxER8SAFuIiIiAf9P5iOXLhkDDzEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In second:\n",
      "Using MAE tau_half: 881062446.551\n",
      "Using MSE tau_half: 881923903.681\n",
      "Real tau_half: 173448000000.000\n",
      "In year:\n",
      "Using MAE tau_half: 27.938\n",
      "Using MSE tau_half: 27.966\n",
      "Real tau_half: 5500.000\n"
     ]
    }
   ],
   "source": [
    "#calculate the predicted tau_half of 20th point\n",
    "y_l1 = predict(x_last, a, b)\n",
    "y_l2 = predict(x_last, a_current, b_current)\n",
    "y_label = y_last\n",
    "W_l1 = 10**y_l1\n",
    "W_l2 = 10**y_l2\n",
    "W_label = 10**y_label\n",
    "tau_half_l1 = np.log(2) / W_l1\n",
    "tau_half_l2 = np.log(2) / W_l2\n",
    "tau_half_label = np.log(2) / W_label\n",
    "#in second\n",
    "print('In second:')\n",
    "print('Using MAE tau_half: {:.3f}\\nUsing MSE tau_half: {:.3f}\\nReal tau_half: {:.3f}'.format(tau_half_l1, tau_half_l2, tau_half_label))\n",
    "#in year\n",
    "conv_rate = 1 / UnitConversion('year')\n",
    "print('In year:')\n",
    "print('Using MAE tau_half: {:.3f}\\nUsing MSE tau_half: {:.3f}\\nReal tau_half: {:.3f}'.format(tau_half_l1 * conv_rate, \n",
    "                                                                                             tau_half_l2 * conv_rate, \n",
    "                                                                                             tau_half_label * conv_rate))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
